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ABSTRACT 

We explore the possibility that the density profiles of elliptical galaxies and cold dark matter (CDM) 
halos found in cosmological simulations can be understood in terms of the same physical process, 
collisionless gravitational collapse. To investigate this, we study a simplified model, the collapse of a 
perfectly cold Plummer sphere. First, we examine an N-body simulation of this model with particles 
constrained to purely radial orbits. This results in a final state characterized by a profile slightly 
steeper than p cx r~ 2 at small radii and behaving as p oc r~ 4 at large radii, which can be understood 
in terms of simple analytic arguments. Next, we repeat our simulation without the restriction of 
radial orbits. This results in a shallower inner density profile, like those found in elliptical galaxies 
and CDM halos. We attribute this change to the radial orbit instability (ROI) and propose a form 
of the distribution function (DF) motivated by a physical picture of collapse. As evidence of the link 
between our model and CDM halos, we find that our collapse simulation has a final state with pseudo- 
phase-space density which scales roughly as p/er 3 oc j- -1-8 , like that observed in CDM halos from 
cosmological simulations (Navarro et al. 2010). The velocity anisotropy profile is also qualitatively 
similar to that found near the centers of these halos. We argue that the discrepancy at large radii 
(where CDM halos scale as p oc r -3 ) is due to the presence of the cosmological background or continued 
infall. This leads us to predict that the outer CDM halo density profile is not "universal," but instead 
depends on cosmological environment (be it an underdense void or overdense region) . 
Keywords: Cosmology: theory — dark matter — Galaxies: elliptical and lenticular, cD 



1. INTRODUCTION 

Over roughly the past 15 years, considerable atten- 
tion has been directed at the possibility that halos in 
cold dark matter (CDM) universes may have so-called 
"universal" density profiles. This idea was motivated by 
studies of cosmological N-body simulations which found 
that the spherically averaged density profiles of CDM ha- 
los could be accurately fit by the same simple fu nction 
for a very large range in mass (jNavarro et al.l ll996) . This 
discovery led to the popular NFW profile parameteriza- 
tion, 



law is given by the Hernquist model (|Hernq uisti fl 9 9 ) 
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where r s is a characteristic length scale. There have 
been numerous attempts to understand the physical 
nature of the NFW profile and other similar param- 
eterizations. The work has mostly focused on an 
origin within the cosmological context, investigating 
the role of large-sc ale structure and cosmo l ogical infall 
and a ccretion (e.g. iLithwick fc Dalalll2011t iDalal et al. 
20101: [Subramanian et al.l l2000t iZukin fc Ber tschingci 



2OlOa0)~ These efforts have generally failed to produce 



a complete physical description of the processes which 
lead to the halo density profiles found in simulations. 

It is well known that elliptical galaxies have surface 
density distributions described by the de Vaucouleurs 
R 1 / 4 law (or its generalization the S ersic law) (e.g. 
iTruiillo et al1l200l IcTaham fc Cuzmanll2003t) . A conve- 
nient analytical approximation to the deprojected i? 1 / 4 
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The NFW and Hernquist profiles are strikingly similar, 
only differing asymptotically at large radii. This suggests 
a link between the density profiles of elliptical galaxies 
and CDM halos. In fact, they can be accurately fit in 
projection by the Sersic law or in three dimensions by a 
related profile, the Einasto model, which has been shown 
to fit CDM ha los more accurately than the NFW pa - 
rameterization (jMerritt et alj|2005t iNavarro et al.ll2010f) . 
This similarity motivates the hypothesis that the same 
physical processes are responsible for the form of the den- 
sity profiles in both CDM halos and elliptical galaxies. 

This leads us to propose that the formation of CDM 
halos can be understood, at least in part, in terms of the 
physics of dissipationless gravitational collapse common 
to elliptical galaxies and CDM halos. We can apply this 
physics to better understand the density profiles observed 
in cosmological simulations. Fo r example, t he a nalytic 
argum ents made 25 years ago by Jaffc (19871) and [White 
(|1987|) imply that in the outer regions of isolated ellipti- 
cal galaxies the density scales as p oc r~ 4 . We suggest, 
that for dark matter halos, any deviation from this is due 
to continued accretion or because at large radii the cos- 
mological background density dominates. Thus, to the 
extent that the structure of halos depends on their in- 
stantaneous accretion state or the local background den- 
sity, the notion that CDM halos have "universal" density 
profiles is not accurate. In a ACDM universe, an isolated 
halo in a very low density environment will stop accret- 
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ing relatively early, and hence will achieve the limiting 
form p oc r -4 . This is con sistent with the results of th e 
cosmological simulations of lDubinski fc CarlbersJ (|199lD , 
where vacuum boundary conditions led to an outer pro- 
file with p oc r~ 4 . We a lso no te that the recent work o f 
iLithwick fc Dalall (|2011[ ) and iVogelsberger et al.1 (|2011[ ) 
demonstrates that for scale-free models, the inner halo 
profile is sensitive to the shape of the collapsing pertur- 
bation. This suggests that the inner profile of halos could 
be sensitive to environment as well. 

In this paper, the first in a series of two, we explore the 
physics of dissipationless gravitational collapse and how 
it applies to the formation of CDM halos. We accomplish 
this with numerical experiments and analytic arguments 
based on a simple toy-model, the collapse of a perfectly 
cold (zero initial kinetic energy) Plummer sphere. In 
Paper 2 (in preparation) we will study the impact of 
cosmological effects, such as continued accretion. This 
will include analyzing cosmological N-body simulations 
to study the expected systematic dependence of the de- 
tailed profile on cosmological environment, be it an un- 
derdense void or an overdense region. 

We begin the study of our collapsing Plummer sphere 
by performing a simulation of purely radial collapse in 
which non-radial motions are suppressed. This results 
in a final equilibrium state with a density profile slightly 
steeper than p oc r~ 2 near the center and p oc r -4 at large 
radii. However, a system composed entirely of particles 
on purely radial orbits will be dynamically unstable once 
non-radial motions are allowed, owing to t he radial orbit 
instability (ROI), an effect fir st noted by iBarnes et al.l 
(p8l and lMerritt fc Aguilarl (p85lh When we remove 
the restriction of purely radial orbits we find a final state 
with a shallower inner profile more similar to that found 
in elliptical galaxies and CDM halos. The possibility that 
the ROI could play a role in structuring CDM halo s has 
been emphasized already by iBellovarv et al.l ([2008) . Be- 
low, we examine this idea in greater depth by postulating 
a physically motivated form for the distribution function 
(DF) of systems which have relaxed through the ROI. As 
evidence of the link between our simple collapse model 
and CDM halos, we find that the final state has a pseudo- 
phase-space density p(r)/a(r) 3 oc r - 1 - 875 j as observed in 
halos in N-body cosmological simulations (jNavarro et al.l 
1201 Of) and der ived in self-similar secondary infall models 
(|Bertschingerlll985fl . 

This paper is structured as follows. In Section 2 we 
present the simulation of our toy-model of gravitational 
collapse with particles restricted to purely radial orbits. 
In the following section, we discuss the outer density pro- 
file of this numerical experiment and how it relates to 
CDM halos and elliptical galaxies. We consider the in- 
ner density profile in Section 4. Here we re-simulate the 
collapse from Section 2 without the restriction of purely 
radial orbits, and discuss the impact of the ROI. In Sec- 
tion 5 we propose a physically motivated DF for the final 
state of systems formed through initially radial gravita- 
tional collapse. We discuss our results in the context of 
previous work in Section 6 and give our main conclusions 
in Section 7. 

Throughout the paper, the data from our numerical 
experiments are presented in the internal un its we use in 
the GADGET N-body code (|Springell [20051) . Distances 



are given in units of lkpc, velocities in units of 1km s _1 , 
and mass in units of 10 10 M Q . Other units (e.g. time, 
energy, density) are given by the relevant combination of 
these three quantities. 

2. A SIMPLE EXAMPLE 

We begin the exploration of our simplified model with 
an N-body simulation of gravitational collapse where par- 
ticle orbits are constrained to be completely radial. The 
simulation was initialized as a perfectly cold (all particles 
with zero velocity) Plummer sphere with scale length of 
100 and total mass of 100 in the code units described 
above. These choices are unimportant, as neither the to- 
tal mass nor the scale length affect the shape of the den- 
sity profile for the final state. The system was evolved in 
time with particles only subjected to the component of 
the gravitational force in the radial direction. This was 
accomplished by modifying the publicly available N-body 
code GADGET (|Springelll2005l ). The simulation was run 
until it reached an equilibrium state, with the density 
profile plotted in Figure [U We have achieved conver- 
gence in the sense that the final state is not sensitive to 
the number of particles included in the simulation. For 
the results shown, we used 10 5 particles which were ini- 
tially distributed by randomly drawing from a probabil- 
ity density function (PDF) proportional to the Plummer 
profile. A gravitational softening length of e g = 0.3 was 
adopted; reducing this value did not affect the results on 
the scales of interest. Note that the specific choice of a 
Plummer sphere does not have a strong impact on the 
final state. We find similar outcomes with other initial 
density profiles (e.g. Gaussian or Hernquist). 

Our simulation results in an equilibrium state slightly 
steeper than p oc r~ 2 at small radii and p oc r -4 at large 
radii. We can show that in general, a density profile 
formed through purely radial collapse cannot be shal- 
lower tha n p oc r~ 2 using fairly simple analytical ar- 
guments (jDalal et al.ll2010f i. In the central regions, the 
time averaged contribution to the total density profile 
from each particle is 
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where dt(r) is the time the particle spends between radii 
r and r+dr. For a density profile shallower than p oc r~ 2 , 
dt/dr is a constant for sufficiently small r. It follows from 
equation ([3]) that central density profiles must always be 
at least as steep as p oc r~ 2 in the purely radial case. 

The large radiu s limit can be explained with the fol- 
lowing argument (jWhitd 119871 : lJaffeHl987l ). In the final 
state of a system formed through gravitational collapse, 
the outer envelope will be populated by particles which 
have been scattered into weakly bound orbits. The size 
of these orbits will be inversely proportional to the en- 
ergy, E. The number of particles between E and E + dE, 
N(E)dE, is expected to be non-zero and continuous near 
E — 0. For the density, this implies 

, , N(E)dE N(GM/r) diGM/r) _ 4 
p[r) = - — T ——(x t. : oc r . (4) 



dr 



Note that this result was derived in the context of el- 
liptical galaxies and has been completely ignored in the 
literature on cosmological CDM halo formation. 
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Figure 1. The density profile (left panel) and logarithmic slope (right panel) of the final state from our simulation of purely radial collapse 
starting from a perfectly cold Plummer sphere (see Section 2 for more details). The solid lines indicate p oc r~ 2 at small radii and p oc r -4 
at large radii. 



Although the above discussion is (intentionally) highly 
idealized, it demonstrates some important consequences. 
First, there is nothing mysterious or pathological about 
density profiles which diverge as r — > 0. Second, the 
large-r behavio r is co nsist ent w i th the analytical argu- 
ments of lWhitd (119871) andlJaffd (11987ft Note that while 
"violent- relaxation" ( Lvnden-Belll 119671 ) plays a role in 
the evolution by redistributing energy among the parti- 
cles, and hence promoting a small fraction to unbound 
orbits, it does not in detail determine the limiting be- 
havior at small and large radii, and hence does not jus- 
tify efforts to understand the nature of the generic den- 
sity profiles based on the rmodynamic arguments (e.g. 
ISpergel fc Hernquistlll992t ). Finally, the lack of sensitiv- 
ity to the number of particles indicates that the outcome 
has been determined by the physics of collisionless dy- 
namics. We argue that the physical processes acting in 
this example are, at least in part, responsible for the ori- 
gin of the density profiles of dark matter halos. While 
the density profile obtained from this example does not 
match the NFW form, we can nevertheless use it as a 
starting point for the more detailed study described in 
the following sections. 

Before moving on, we note that the final state in this 
example can be approximated by the density profile of 
the Jaffe model, 
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However unlike the isotropic Jaffe model, the velocities 
are purely radial and therefore are dynamically unstable 
when particles are permitted to have non-radial velocities 
(as we explore below). Note that the DF for a Jaffe model 
density profile constrained to purely radial orbits can be 
expressed analytically. We derive this expression in the 
Appendix. 

3. THE OUTER DENSITY PROFILE 

Acc ordin g to t he ar guments put forward by I White! 
(|1987ft and I Jaff e (1987), an isolated collisionless system 
evolving into equilibrium through collapse should even- 
tually achieve the limiting form p oc r -4 at large radii. 
This expected outcome is in agreement with the simple 
numerical experiment presented in Section 2. Note that 
the analytic argument is insensitive to orbital anisotropy, 



and therefore is unaffected by the imposed constraint of 
purely radial orbits. We point out that this large radius 
behavior should also occur i n dark matter halos once cos - 
mological infall terminates. Dubinski fc Carlbergl ([1991ft 
originally demonstrated this using cosmological simula- 
tions in which they employed vacuum boundary condi- 
tions around individual objects which naturally led to a 
cessation of late time accretion. 

We suggest that the NFW profile at large radii, p oc 
r -3 , does not directly reflect the physical processes that 
form CDM halos. We find in our collapse simulations 
(both the purely radial case discussed above and the 
generalized case discussed below) that the density pro- 
file reaches p oc r~ 4 roughly two orders of magnitude 
beyond the radius where cfln(p) j d\n(r) ~ —2. For most 
CDM halos in the cosmological context, the background 
density will dominate at such large radii. Thus, CDM 
halos may not extend far enough before falling below the 
background density to achieve p oc r~ 4 . Additionally, 
continued infall in high density regions could alter the 
outer profile. This finding has important implications 
for understanding the formation of halos in a cosmolog- 
ical context. Low-density regions will behave as open 
universes, implying that accretion there will halt ear- 
lier than in higher density environments. Some regions 
may be underdense enough to permit halos to reach the 
asymptotic form p oc r~ 4 . Therefore, it is expected that 
the outer profiles of halos may exhibit a variety of be- 
haviors, depending on environment and redshift, so in 
detail the outer profiles of dark matter halos will not be 
characterizing by a universal form. This will be explored 
in detail in paper 2 (in preparation). 

4. THE INNER DENSITY PROFILE AND THE RADIAL 
ORBIT INSTABILITY 

To investigate the small-radii properties of equilibrium 
systems formed through gravitational collapse, we gen- 
eralize the numerical experiment described in Section 2. 
We perform the same N-body simulation of radial col- 
lapse, but relax the restriction that the particles move 
on purely radial orbits. Here we use 2 x 10 6 particles 
and a softening length of e g = 0.1. We find that we have 
achieved convergence in the sense that increasing particle 
number or decreasing softening length do not impact our 
results on the scales of interest. Since our results are not 
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sensitive to particle number, we conclude that 2-body 
relaxation does not impact our results and the system 
is effectively collisionless. As in the purely radial sim- 
ulation, different initial density profiles result in similar 
final equilibrium states. The results of this simulation are 
presented in Figure [21 Note that the simulation was run 
roughly five times the dynamical time at the maximum 
radius shown. While we may not have reached perfect 
equilibrium at this radii, we do not expect substantial 
changes at longer times. 

Initially the system has mainly radial orbits (as in 
the first numerical experiment). However, in response 
to small non-radial perturbations, purely radial orbits 
are dynamically unstable. The system achieves a differ- 
ent equilibrium in which the orbits in the inner regions 
are isotropic, but those in the outer regions, where the 
system is not strongly self-gravitating, are still mainly 
radial. Roughly, the transition occurs at the half-mass 
radius of the new configuration. This is shown in Fig- 
ure [2 where we plot the radial profile of the anisotropy 
parameter, 



where a r (r) 2 and <J t (r) 2 are, respectively, the radial and 
the tangential variance in the velocity as a function of 
radius ((3 — corresponds to isotropic and = 1 to ra- 
dial). In the final state, the system has an inner density 
profile shallower than r~ 2 , but the outer regions are un- 
changed compared to the purely radial case, still having 
p oc r~ 4 . The particles in this outer region (well out- 
side the half-mass radius) respond to the gravitational 
potential dictated by the inner region, and their orbits 
are therefore stable. The inner core (inside the half mass 
radius), on the other hand, is self-gravitating and sub- 
ject to the ROI. The resulting central density cusp is 
now more similar to those which characterize CDM ha- 
los and the deprojected i? 1 / 4 law, supporting the notion 
that CDM halos and elliptical galaxies are shaped by 
similar processes. 

As further evidence of the relationship between the 
physics of gravitational collapse in this example and 
CDM halos we observe that the pseudo-phase-space den- 
sity of the final state scales roughly as p/a 3 cx r~ 1875 , 
matching that in ha los from cosmological simulations 
(jNavarro et all 120100 . This is plotted in Figure O 
We also find that the velocity anisotropy profile scales 
roughly as /3(r) oc ln(r) near the center, qualitatively 
similar to that seen in halos from cosmological N-body 
simulations. 

5. A DISTRIBUTION FUNCTION FOR THE RELAXED 
SYSTEM 

The numerical experiments and analytic arguments 
presented in the preceding sections have led us to propose 
a simple scenario for how collisionless systems are formed 
through gravitational collapse. The outer profile is de- 
termined by the redistribution of energy among particles 
through violent-relaxation, which promotes a small frac- 
tion of the mass onto unbound orbits. If no further infall 
occurs, the density profile will asymptote t o p oc r -4 as 
predicted bv I White! (fl987h and[Ja|l <fl987h . If infall has 
not halted, the system will not be in equilibrium and a 
different density profile may be measured at large radii. 



For a nearly radial infall, the ROI will operate causing 
the system to evolve such that the inner region is charac- 
terized by orbit isotropy and a density profile shallower 
than p oc r" 2 (not permitted in the purely radial case). 
If the collapsing particles initially have significant non- 
radial motions, the final profile will be even shallower. 

To understand the physical nature of systems formed 
through this scenario, we have performed an analysis of 
the phase-space distribution of the final state in the nu- 
merical example of Section 4. Note that because the 
system is anisotropic, the DF cannot be expressed as a 
function of energy alo ne. We propose the following form 
of the DF (guided bv lGerhardl (fl99l 'l. 

f(E,L) = f(E,x) = g(E)h(x,E), (7) 

where x = L/L C (E) is the angular momentum divided by 
the angular momentum of a particle on a circular orbit 
with energy E. The values of x range between and 1. 
The functions g and h are given by 

g{E) = A C M~P E E)(~Ef 1 (8) 

and 

h{x,E)=e^>{-x 2 /q{E) 2 ), (9) 

where A, and c are free parameters. We expect q(E) 
to be a function which increases with \E\ and parame- 
terize it with a broken power law as explained below. 

This model's main features can be explained in terms 
of the collapse scenario described above. The function 
h(x, E) determines how anisotropic orbits of a given en- 
ergy are, with a delta function in x corresponding to 
purely radial and a constant function corresponding to 
isotropic. We expect that q(E) gets larger for increasing 
\E\ such that h(x, E) transitions from a delta function 
for the least tightly bound particles to a constant for the 
most tightly bound particles. The most tightly bound 
particles are those which make up the inner profile and 
have been isotropized through the ROI, while the more 
weakly bound particles are unaffected by this instability 
and are found mostly in the outskirts of the system on 
radial orbits. The energy dependence of the DF in equa- 
tion (jHJ is a Maxwellian distribution for large \E\, but 
modified by the factor (—E) c for the least tightly bound 
particles. We can interpret this as a Gaussian distribu- 
tion (representing a favorable statistical state), modified 
at high energies due to particles being scattered into un- 
bound orbits. The values of the parameters in our model 
can also be related to the bulk properties of the system. 
Our expectation is that the parameter l//3 e is of the or- 
der of the kinetic energy per particle. The transition 
in q(E) which changes h(x,E) from a delta function to 
a constant should happen roughly at the mean energy 
of particles found at the half-mass radius of the system. 
This corresponds to the notion that the ROI isotropizes 
roughly the inner half of the system and leaves the outer 
envelope, which is not strongly self-gravitating, mostly 
radial. 

We compare the DF measured directly from the sim- 
ulation with a parameterization of this model. Values 
of c — 0.9275 and f3 E = 4 x 10~ 5 are found to match 
the simulation well (see Figure |4]). As expected, the 
value of 0e also corresponds roughly to the reciprocal 
of the mean kinetic energy, v 2 /2 = 10 4 . For q(E) we 
use a broken power law of the form q(E) — a(—E) h . 
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Figure 2. The final state of the N-body simulation of collapse from a cold (zero velocities) Plummer sphere described in Section 4. We 
plot radial profiles of the density (upper left), cumulative mass interior to a radius r (upper right), anisotropy parameter (lower left), and 
logarithmic slope (lower right). The solid line indicates p oc r -4 and the dashed line marks the half-mass radius. 

tion and this parameterization of the model in Figure 2J 
Clearly, there is very good agreement. 

In addition to fitting the phase-space data from the 
numerical experiment we also take this DF parameter- 
ization and self-consistently compute the density and 
anisotropy profiles of the system. This was accomplished 
iteratively with the following expression: 
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Figure 3. The pscudo-phase-space density, p/o" 3 , of the final state 
of the collapse simulation described in Section 4 (points) . For ref- 



erence we plot p/cr 



which is found in CDM halos from 



cosmological N-body simulations. 

For ln(-.E) < 9.55, a = 0.00012 and b = 0.87, while 
for 9.55 < ln(-E) < 11.9, a = 0.07 and b = 0.21. When 
ln(-_B) > 11.9, we assume a = 4x 10~ n and 6 = 2. This 
was determined by fitting for q(E) separately in each en- 
ergy bin and then using the power law as a simple way 
to parameterize the values. This procedure provides a 
reasonably good match to the simulation data and has 
the limiting behavior we expect for high and low E. We 
compare contours of the DF, f(E,L), from the simula- 
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(r) 



d 3 vf(E = v 2 /2 + <f>(r),L = V t r), (10) 



where the gravitational potential $(r) is computed using 
the density from the previous step, p n . At each step the 
normalization of the DF is set to fix the innermost value 
of p(r) to a constant value. We use the density profile 
from the simulation as p(r)i and iterate until the pro- 
file has converged to values where p{r) n = p(r) n +i. In 
Figure [5j we show the results of this procedure and find 
good agreement with the density and anisotropy profiles 
measured from the simulation. Note that we have not 
searched the parameter space for the best match between 
profiles. We simply used the same parameters described 
above. A better match would be obtained if we searched 
explicitly for the values of c, Pe, and the q(E) power law 
parameters which minimize the density profile discrep- 
ancy. More flexible functions of q(E) could also improve 
the fit. 

That a relatively simple form of the DF matches the re- 
sults of our numerical experiments so well justifies the no- 
tion that there is favored outcome through gravitational 



6 



collapse. A similar outcome should apply to all systems 
formed through radial dissipationless gravitational col- 
lapse, including CDM halos and elliptical galaxies. 

6. COMPARISON WITH PREVIOUS WORK 

There has been much work attempting to understand 
the physical nature of the nearly universal density pro- 
files of CDM halos found in cosmological N-body simula- 
tions. Next we summarize much of this work and explain 
how it relates to the present paper. 

One popular way of analytically modeling CDM halo 
formation is through so-calle d secondary infall (SI) mod- 
els. This work began with iGunn fc Gottl (|1972l i. who 
studied how mass shells accrete onto a collapsed ob- 
ject. Later, by im posing s elf-similar and radia l orbits 
iBertschingen (p85T ) and iFillmore fc Goldreichl (fl98l 
were able to derive the inner profile of a halo as a func- 
tion of the density profile of the initial perturbation. 
More recently, there have been a number of improve- 
ments to these models which i nclude angular momentum 
(|Zukin fc B crtschingcr 2010 aH) and generalize to fully 
3D solutions ( Lithwick fc Dalall2011f), incl uding numer- 
ical simulations (jVogelsberger et al.l feplll ). Note that 
much like our example in Section 2, the purely radial SI 
models produce halos which have steeper central profiles 
than those observed in cosmological simulations. Models 
which include non-radial velocities are able to produce 
NFW or Einasto-like profiles. However, they contain a 
number of non-physical assumptions to make the prob- 
lem tractable, such as spherical symmetry or assuming 
that all particles in a given mass shell have the same ve- 
locity. In this paper, we argue that much of the physics 
which determines the density profiles of CDM halos is 
unrelated to cosmology. The details related to cosmolog- 
ical expansion and accretion included in SI modeling may 
obscure the important physics. Thus, our work compli- 
ments the SI approach by identifying what aspects of the 
equilibrium state do not depend on cosmology. 

Other work focusing on the cosmological context has 
investigated how the hier archical assembly of halos af- 
fects the density profile. iHuss et all (j!999h found that 
a variety of different initial conditions, including sim- 
ple monolithic c ollaps e produce similar density profiles. 
iWang fc White] ([2009D find that cosmological hot dark 
matter (HDM) simulations produce similar profiles to 
CDM halos, implying that the mergers due to increased 
substructure in a CDM cosmology are not important to 
structuring density profiles. This is a lso su pported by 
the simulations of iKazantzidis et ail (|2006f) , who find 
that the shape of density profiles are preserved during 
major merg ers of halo s in N - body simulations, and ear- 
lier wo rk bv lHernqu ist (1992, 19931): [Barnes fc Hernquistl 
(1996) who showed that density profiles similar to those 
described here arise naturally in mergers of individual 
galaxies. All of this work supports the idea presented in 
this paper, that the shape of CDM halos is due to physi- 
cal principles common to collapsing systems rather than 
depending on the details of initial conditions or cosmol- 
ogy- 
It has also been observed in cosmological simula- 
tions that the outer portions of CDM halo density pro- 
files depend on the cosmological model. For example 
the outskirts of halos simulated with vacuum boundary 
conditions have profiles steeper than the NFW profile 



(jDubinski fc CarlberdH99l . In the future of a ACDM 
universe, halos are steeper than the NFW parameteriza- 
tion at large radius and are trunc ated due to t he ac - 
celeration caused by dark energy (|Busha et al.ll2005h . 
We note that the analytic prediction of p oc r -4 breaks 
down at large radii in this case since the binding energy 
of particles is affected by dark energy. That the outer 
profile is affected by cosmology supports our claim that 
the p oc r~ 3 asymptotic limit in the NFW profile is a 
consequence of the cosmological background density and 
continued accretion and not due to the primary physics 
which determines the profile. 

We argue that the ROI is a key physical process which 
determines the shape of CDM halo density profiles. This 
has been proposed by a number of other authors and 
the R OI has been explored both within N-body s imula- 
tions (iBellovarv et all 120081: iMacMillan et al.ll2006D and 
semi- analytic models ([Barnes et al.l 120051 ) . We extend 
these studies by proposing a form of the DF motivated 
by a physical scenario of gravitational collapse including 
the ROI. A different form of the DF to describe objects 
formed through gravita tion al collapse was put forward by 
iTrenti fc Bertinl (|2005l ) and lTrenti etTal] (|2005l ). This DF 
was determined by maximizing entropy subject to a con- 
straint involving an integral of the angular momentum 
and the energy. This model generates density and veloc- 
ity anisotropy profiles similar to those found in N-body 
simulations, but has a feature which does not match in 
the detailed f o rm of f(E,L). Another relevant paper is 
iMerritt et al.l (|2005f ). where they point out, as we have, 
that elliptical galaxy and CDM halo profiles can both 
be fit by similar parameterizations (namely the Einasto 
profile or Sersic model). In this paper we have attempted 
to better explain this connection by studying a new form 
of the DF which better matches simulations. 

7. CONCLUSIONS 

The similarity between the density profiles of CDM 
halos and elliptical galaxies suggests that their forma- 
tion can be understood in terms of the same physical 
processes. This leads us to propose that the features 
of CDM halo density profiles result, to a large extent, 
from the physics of gravitational collapse of collisionless 
systems. In this paper, the first in a series, we studied 
these dynamics by exploring a toy-model, the collapse of 
a perfectly cold Plummer sphere. 

Using analytic arguments and two numerical experi- 
ments, we demonstrated that the final equilibrium state 
of our Plummer model collapse can be understood in 
fairly simple terms. In the first experiment, we ran an N- 
body simulation of collapse with the particles constrained 
to radial orbits. In the final state, we have found an in- 
ner density profile which is slightly steeper than p cx r~ 2 
and an outer profile which scales as p oc r~ 4 . Sim- 
ple analytic arguments explain these asymptotic profiles. 
In the second experiment, we relaxed the radial orbit 
constraint. Here we have found the same outer equi- 
librium profile, but a shallower inner profile, which we 
attribute to isotropization introduced by the ROI in the 
self-gravitating core of the halo. 

We are led to a relatively simple scenario for the for- 
mation of objects through radial gravitational collapse. 
Initially violent-relaxation scatters some particles to un- 
bound orbits and the remaining weakly bound particles 
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Figure 4. A comparison of the DF from the N-body simulation of gravitational collapse described in Section 4 and the model DF described 
in Section 5. We plot contours of the DF, log(/(_E, L)), from the simulation (upper left panel) and the parameterized model (upper right 
panel). In the lower left panel we plot g(E), given by equation (8), for the simulation data (blue points) and the model (dashed line). In 
the lower right panel we plot h(x, E = —3.8 X 10 4 ), given by equation JS}. The values for the model parameters are listed in the main text. 





Figure 5. A comparison between the numerical simulation described in Section 4 and the model DF described in Section 5. We plot the 
density profile (left panel) and the velocity anisotropy (right panel) from the simulation (blue points) and the values implied by the model 
DF (dashed lines). The values for the model parameters used are listed in the main text. 



form the outer d ensity profile , whi ch scale s as r 4 due to 
the argument of Uaffel (fl987[ ) and iWhitd (|1987l) . These 
particles are not strongly self-gravitating and thus re- 
main on mostly radial orbits. The inner profile is strongly 
self-gravitating and becomes isotropic through the onset 
of ROI resulting in a profile shallower than p oc r~ 2 , like 
those observed in elliptical galaxies and CDM halos. 

We present a form of the DF which matches the final 
state of our collapse simulation. The main features of this 
model can be related to our physical picture of collapse; 



weakly bound (large radii) particles are on radial orbits 
and tightly bound (small radii) particles have isotropic 
velocities. This is captured with a term that describes 
how circular orbits are as a function of energy, h(x,E). 
These simple considerations are combined with a roughly 
Maxwellian energy dependence term, which is a natural 
expectation based on statistical arguments. Our model 
self-consistently reproduces the density and anisotropy 
profiles obtained from our simulations. The parameters 
in the model can also be related to the bulk properties 



of the system. 

The simple parameterization of the DF motivates a 
universal form for density profiles created through this 
process. However, non-radial velocities in the initial con- 
ditions of collapse may alter the inner profile, making 
it shallower than we find in our numerical experiments. 
Continued cosmological infall and accretion could also 
have some effect. These factors may make the inner pro- 
files of CDM halos sensitive to cosmological parameters 
in detail. Since the outer profile of an object created 
through gravitational collapse is predicted to scale as 
p oc r~ 4 , the shallower profile in CDM halos must be 
due to ongoing cosmological accretion and the presence 
of the cosmological background. 

As evidence for the connection between our toy-model 
and CDM halos we find that they both have the same 
power law for the pseudo-phase-space density, p/er 3 oc 
r -1,875 . We also note that at small radii, CDM halos 
have an anisotr opy profile, 0(r), wh ich is similar to our 
collapse model (|Navarro et al.|[20~10h . 

The physical picture presented in this paper provides 



predictions for CDM halos which can readily be tested. 
First, there is no "universal" form of the CDM halos in 
the Universe. Instead, the outer halo is altered by accre- 
tion and therefore will depend on background density. In 
principle, a halo in a sufficiently under-dense region will 
behave as if it is in an open-universe, causing accretion 
to halt at earlier epochs allowing the halo to equilibrate 
to the p oc r~ 4 form. Another prediction is that the final 
state of halos may have DFs which resemble those from 
our collapse simulation. We will test these predictions in 
paper 2 (in preparation) and further investigate how the 
physics of dissipationless collapse and effects from the 
cosmological environment combine to produce the CDM 
halo density profiles seen in N-body simulations. 
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APPENDIX 

DISTRIBUTION FUNCTION OF THE RADIAL JAFFE PROFILE 
The distribution function (DF) for a system of purely radial orbits can be written as 

f(x,v)=g(£)S D (L 2 ), 



(Al) 



where £ = —E and S D is a Dirac delta function. We derive g(£ ) for a Jaffe density profile with scale length a and 
density parameter po, 

P(r) = , , ^ , (A2) 



(r/a) 2 (l + rjaf 



We begin with 



p(r) = J d 3 vf(x, v) = J dv, 



where v r is the radial velocity. We then perform a change of variables to obtain 

"* d£g(£) 



p{r) 



o x/* - £' 



(A3) 



(A4) 



where we have defined p (r) = \/2p(r)r 2 and = — $. Using the Abel integral equation (see Appendix B in 
Bin nev fc Tremai ne (2008)), we can solve for g(£), 



9(£) = z 



dp p(* = 0) 



o V£ - # d * V£ 



Since ^/ is a monotonic function of r we are able to use p as a function of ^P, 

2 



2„-2* / „* 



1 



where ^ = ^//(47rpoa 2 G). Evaluating the integral in equation (| A5|) we obtain, 

2a [po 



7T 3 / 2 V G 



V2D 



I £ 
4ira 2 p G 



D 



£ 



2na 2 poG 



(A5) 



(A6) 



(A7) 



where D(x) is the Dawson function. 



REFERENCES 



Barnes, E. I., Williams, L. L. R., Babul, A., & Dalcanton, J. J. 
2005, ApJ, 634, 775 



9 



Barnes, J., Hut, P., & Goodman, J. 1986, ApJ, 300, 112 
Barnes, J. E., & Hernquist, L. 1996, ApJ, 471, 115 
Bcllovary, J. M., Dalcanton, J. J., Babul, A., et al. 2008, ApJ, 
685, 739 

Bertschinger, E. 1985, ApJS, 58, 39 

Binney, J., & Tremainc, S. 2008, Galactic Dynamics: Second 

Edition (Princeton University Press) 
Busha, M. T., Evrard, A. E., Adams, F. C, & Wcchsler, R. H. 

2005, MNRAS, 363, Lll 
Dalai, N., Lithwick, Y., & Kuhlen, M. 2010, ArXiv e-prints 
Dubinski, J., & Carlberg, R. G. 1991, ApJ, 378, 496 
Fillmore, J. A., & Goldrcich, P. 1984, ApJ, 281, 1 
Gerhard, O. E. 1991, MNRAS, 250, 812 

Graham, A. W., & Guzman, R. 2003, The Astronomical Journal, 
125, 2936 

Gunn, J. E., & Gott, III, J. R. 1972, ApJ, 176, 1 
Hernquist, L. 1990, ApJ, 356, 359 
— . 1992, ApJ, 400, 460 
— . 1993, ApJ, 409, 548 

Huss, A., Jain, B., & Steinmetz, M. 1999, ApJ, 517, 64 
Jaffe, W. 1987, in IAU Symposium, Vol. 127, Structure and 

Dynamics of Elliptical Galaxies, cd. P. T. de Zeeuw & S. D. 

Tremainc, 511 

Kazantzidis, S., Zentner, A. R., & Kravtsov, A. V. 2006, ApJ, 
641, 647 

Lithwick, Y., & Dalai, N. 2011, ApJ, 734, 100 
Lynden-Bell, D. 1967, MNRAS, 136, 101 



MacMillan, J. D., Widrow, L. M., & Henriksen, R. N. 2006, ApJ, 
653, 43 

Merritt, D., & Aguilar, L. A. 1985, MNRAS, 217, 787 
Merritt, D., Navarro, J. F., Ludlow, A., & Jenkins, A. 2005, ApJ, 
624, L85 

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 
563 

Navarro, J. F., Ludlow, A., Springel, V., et al. 2010, MNRAS, 
402, 21 

Spergel, D. N., & Hernquist, L. 1992, The Astrophysical Journal 

Letters, 397, L75 
Springel, V. 2005, MNRAS, 364, 1105 

Subramanian, K., Ccn, R., & Ostriker, J. P. 2000, ApJ, 538, 528 
Trenti, M., & Bertin, G. 2005, A&A, 429, 161 
Trenti, M., Bertin, G., & van Albada, T. S. 2005, A&A, 433, 57 
Trujillo, I., Erwin, P., Asensio Ramos, A., & Graham, A. W. 

2004, The Astronomical Journal, 127, 1917 
Vogelsberger, M., Mohayaee, R., & White, S. D. M. 2011, 

MNRAS, 414, 3044 
Wang, J., & White, S. D. M. 2009, MNRAS, 396, 709 
White, S. D. M. 1987, in IAU Symposium, Vol. 127, Structure 

and Dynamics of Elliptical Galaxies, ed. P. T. de Zeeuw & 

S. D. Tremaine, 339-349 
Zukin, P., & Bertschinger, E. 2010a, Phys. Rev. D, 82, 104044 
— . 2010b, Phys. Rev. D, 82, 104045 



